StatBinscatter <- ggplot2::ggproto(
  "StatBinscatter", 
  Stat,
  compute_group = function(data, scales, bins = 10) {
    bins     <- min(floor(nrow(data)/10), bins)
    x_bin    <- ggplot2::cut_number(data$x + 1e-12*runif(nrow(data)), bins)
    x_means  <- stats::ave(data$x, x_bin, FUN = mean)
    y_means  <- stats::ave(data$y, x_bin, FUN = mean)
    y_se     <- stats::ave(data$y, x_bin, FUN = sd)
    y_obs    <- stats::ave(data$y, x_bin, FUN = length)
    result   <- data.frame(x    = x_means, 
                           y    = y_means, 
                           ymax = y_means + 1.96*y_se/sqrt(y_obs),
                           ymin = y_means - 1.96*y_se/sqrt(y_obs))
    result   <- unique(result)
    return(result)
  },
  required_aes = c("x", "y")
)

stat_binscatter <- function(mapping = NULL, data = NULL, geom = "point",
                            position = "identity", na.rm = FALSE, show.legend = NA, 
                            inherit.aes = TRUE, ...) {
  ggplot2::layer(
    stat = StatBinscatter, data = data, mapping = mapping, geom = geom, 
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
    )
}

# These codes are cited from Zeng et al (2023)
# Zeng, Y., Jiang, J., Li, J., & Göbel, C. (2023). The rise of grassroots civil society under one-party rule: The case of China’s homeowner associations. World Politics, 75(3), 608–646. https://doi.org/10.1353/wp.2023.a900714
